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COMPUTER  SIMULATION  STUDIES  OF  THE  TEARING  MODE 
INSTABILITY  IN  A  FIELD-REVERSED  ION  LAYER 

I.  INTRODUCTION 

The  advent  of  pulsed,  high  power  ion  sources"''  during 

the  last  several  years  has  stimulated  renewed  interest  in 

the  formation  of  field-reversed  configurations  (FRC)  with 

ion  rings  or  layers.  Similarly,  the  availability  of  intense, 

long  duration  neutral  beams  has  focused  attention  on  the 

possibility  of  forming  field-reversed  configurations  in  an 

2 

ordinary  magnetic  mirror  by  either  ion  layers  or  plasma 

currents.3  Presently,  thermonuclear  reactors  based  on  such 

field- reversed  configurations  look  very  attractive  because 

their  anticipated  high  power  density  might  lead  to  relative 

4 

low  cost  and  modest  size  systems. 

The  usefulness  of  ion  ring  or  layer  generated  FRC  rests 
very  heavily  on  their  stability.  As  a  result  of  the  complexity 
of  the  problem,  no  comprehensive  treatment  of  the  stability  of 
such  states  is  presently  available.  However,  all  the  results 
obtained  so  far  are  very  encouraging.  Using  the  energy  prin¬ 
ciple  Sudan  and  Rosenbluth5  have  shown  with  a  hybrid  model 
that  a  field-reversed  ion  layer  immersed  in  a  background 
plasma  is  stable  against  MHD  kink  modes  when  its  aspect  ratio 
(=  major  radius/minor  radius)  is  near  unity.  Lovelace^  had 
arrived  earlier  at  a  similar  conclusion.  Uhm  and  Davidson^ 
have  treated  the  stability  of  a  field-reversed  layer  for 
frequencies  near  multiples  of  the  mean  rotational  frequency 
Manuscript  submitted  June  30,  1980 
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of  the  layer.  They  have  found  that  these  instabilities, 
which  are  of  primary  concern  can  be  easily  stabilized  by 
introducing  transverse  temperature  in  the  rotating  ions  of 

g 

the  layer.  Finn  and  Sudan  have  investigated  the  effect  of 
particle  resonances  on  low-frequency  magnetohydrodynamic 
modes  in  a  field  reversed,  long  ion  layer. 

The  tearing  modes,  i.e.,  low  frequency  instabilities 

that  break  up  a  very  long  layer  into  several  rings  have  been 

studied  by  Marx9  and  recently  by  Uhm  and  Davidson.10  Fowler11 

has  found  that  finite  length  layers  are  stable  to  tearing 

modes  as  a  consequence  of  axial  kinetic  pressure.  Using  a 

hybrid  model,  in  which  the  ion  layer  is  treated  kinetically 

and  the  background  plasma  as  a  fluid,  Uhm  and  Davidson  extended 

the  calculation  of  Fowler  to  infinitely  long  layers  and  have 

shown  that  when  the  axial  velocity  spread  Av  exceeds  a 

critical  value  Av  . the  layer  does  not  break  into  rings, 

crxt 

i.e.,  is  stable  to  the  tearing  mode. 

In  this  paper,  we  report  results  from  the  numerical 
simulation  of  an  infinitely  long,  space-charge  neutral  ion 
layer.  The  numerical  results  give  a  ^vcrj_t  that  is  in 
excellent  agreement  with  the  theory.  However,  the  growth 
rates  predicted  by  the  theory  are  always  greater  by  a  factor 
of  two  than  those  found  from  the  code.  An  interesting  result 
of  the  simulation  is  that  although  close-fitted  metal  walls 
always  improve  the  stability  of  an  unstable  layer,  there 
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are  conditions  under  which  a  conducting  shell  that  is  not  very 
near  the  outer  edge  of  the  layer  can  have  an  adverse  effect  on 
the  stability  of  the  layer. 

In  Section  II  of  this  paper  we  give  a  brief  description 
of  the  particle-in-cell  code  used  for  these  studies.  The 
initiation  of  the  rigid  rotor  equilibrium  with  finite  tem¬ 
perature  is  also  described  as  are  the  basic  diagnostics  for 
determining  the  energy  of  the  magnetic  field  spectrum.  The 
linear  theory  for  the  onset  of  the  tearing  mode  is  reviewed 
in  Section  III.  It  is  clear  from  the  interpretation  of  the 
threshold  condition  that  there  exist  regimes  where  the  presence 
of  a  conducting  wall  can  have  a  destabilizing  influence  on  the 
tearing  mode.  In  Section  IV  we  present  the  results  of  the 
simulation  which  confirm  the  existence  of  the  threshold  con¬ 
dition  as  determined  by  the  thermal  spread  in  the  parallel 
velocity  Vz  and  the  geometric  position  of  the  beam. 

The  effect  of  wall  destabilization  is  also  demonstrated 
as  discussed  in  the  conclusion  presented  in  Section  V. 

II.  DESCRIPTION  OF  THE  CODE 

a)  Field  Solver 

The  dynamics  of  reversed  magnetic  field  ion  layers  are 

studied  in  a  cylindrical,  azimuthally  symmetric  (r,z)  geom- 

12 

etry  using  a  2D-3V  particle  simulation  code.  The  annular 
proton  layers  are  space-charge  neutralized  by  electrons  which 
follow  the  radial  and  axial  motion  of  ions  thus  cancelling 
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the  ion  current  in  these  directions.  The  electrons,  however, 
do  not  cancel  the  azimuthal  ion  current.  Thus,  Maxwell’s 
equations  reduce  to 


Ve 
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where  A0  (r,z,t)  is  the  magnetic  vector  potential,  J0  (r,z,t) 
is  the  self  consistant  ion  current  and  is  the  permeability 
of  free  space. 

For  this  azimuthally  symmetric  system  Bz  and  Br  are 
obtained  by  differentiating  A0 .  The  canonical  angular  momentum 
P,  =  r(m.V„  +  qA„)  is  a  constant  of  the  motion  and  is  used  to 
find  V0  of  each  ion. 

The  magnetic  vector  potential  A0  is  solved  on  a  two 
dimensional  (r,z)  grid  (typically  32  x  32  but  expandable  to 
32  x  128)  by  Fourier  analysis  in  the  axial  direction  and  by 
Gaussian  elimination  in  the  radial  direction.  The  system  is 
periodic  in  z  and  is  enclosed  radially  by  a  perfectly  con¬ 
ducting  wall.  The  standard  five  point  difference  operator 
is  used  to  represent  the  differential  operator  V2 . 

The  equations  of  motion  are  solved  in  a  leapfrog  scheme. 
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where  AT  is  the  unit  timestep  and  T  is  the  time.  The 

fields  B  ,  Br  and  vector  potential  A0  are  found  at  the 
i  ri  i 

particle  positions  by  volume  averaging  the  fields  at  the 
four  nearest  cells. 

The  field  solver  and  particle  pusher  are  not  naturally 
time  centered.  The  lack  of  centering  can  and  does  cause  the 
solution  of  the  field  solver  to  become  numerically  unstable. 


To  obtain  a  numerically  stable  solution  we  rewrite  the 


differential  equation  for  AQ  as 
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where  ATV.  “  =  AQ  ^(T  +  3/2)  At)  is  the  new  solution, 


.T+3/2 

VK+1 


A-3/2  _  ^(T  +  1/2)  At^  is  the  solutibn  from  the  previous 

time  step, 

V 

2  A, 


tT+3/2  _ 

JK  “  J9 


((T+3/2)  At)  =  ne  2  (jii  - 


m  k,  i) 


a  is  a  suitably  chosen  constant 

PQ  is  the  canonical  angled  momentum  of  the  ith 
particle , 

r^  is  the  radial  position  of  i  particle, 

T+3/2  th 

Ak  i  =  ^  evaluated  at  the  position  of  the  i 

particle  and  K  is  the  iteration  number. 

The  equation  can  be  written  as 


AK+1  ~  2ctAK+l 


=  2Jr  -  V2Ak  -  2aAK 
=  2Jk  -  (?2 Ak  -  2aAK)  -  4aAK 


with 

SK+1  2  JK  ~  SK  4aAK* 

Thus  for  the  iteration  one  needs  to  retain  only  SR  and  AR 
from  previous  iteration  and  evaluate  J„  after  each  iteration. 
A  discussion  of  the  properties  and  stability  of  this  scheme 
will  be  given  in  another  paper. 

b)  Initialization  of  Particles 
In  the  numerical  experiments  the  macroparticles 
representing  the  protons  are  loaded  to  satisfy  the  equilibrium 
state  described  by  the  distribution  function 


f (H  -  wPa,v_)  = 

v  <u 

X  6  (H  -  0)PA  -  K)  9 [  ( v  +  Av  )  (Av  -  v  )], 
Q  0  z  z  z  z 

where  H  is  the  hamiltonia  perpendicular  to  the  external 
magnetic  field,  XQ,  o  and  K  are  constants,  P0  is  the 
canonical  angular  momentum  and 

0  (x)  =  j1,  X  -  0 

!o,  x  <  o. 
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The  ion  layer  that  is  described  by  the  distribution  function 

1 3 

of  Eq.  (1)  is  in  a  rigid  rotor  equilibrium  with  zero  average 
axial  velocity  and  with  a  total  axial  velocity  spread  of  2  Avz . 

The  evolution  of  the  system  is  determined  by  following 
the  trajectories  of  6000  to  24000  macroparticles  in  the  exter¬ 
nally  applied  magnetic  and  self-consistent  magnetic  fields 
produced  by  the  particles. 

The  distribution  function  of  Eq.  (1)  requires  that  at 
t  =  0  the  following  constraints  be  satisfied: 

<  v  >  =  Id3  r  v  f  =  0 

z  J  z 

<  VR>  =  J d 3  r  vR  f  =  0 

and 

<  vg(r)  >  =  j dzrvgf  =  uir.  (4) 

To  insure  that  the  conditions  of  Eqs .  (2)  -  (4)  are  satisfied, 
an  octet  of  particles  is  loaded  at  each  initializing  position. 

A  random  value  of  v  is  selected  between  0  and  iv  .  Four  of 

2  Z 

the  particles  are  given  the  velocity  Avz  and  the  other  four 
are  given  -  Avz.  The  radial  velocity  is  assigned  in  a  similar 
manner.  The  final  component  Vg  is  given  by  adding  and  sub¬ 
tracting  a  random  velocity  from  oar.  The  result  is  that  each 
of  the  eight  particles  in  an  octet  is  given  in  a  different 
position  in  velocity  space.  The  positions  of  the  particles 
in  the  octet  are  then  randomized  within  a  specified  radial  region. 


(2) 

(3) 
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c)  Diagnostics 

One  of  the  main  diagnostics  used  to  determine  whether 


the  system  is  stable  or  unstable  is  the  evaluation  of  the 
equipartion  of  the  magnetic  energy  of  the  system  in  the  axial 
Fourier  modes.  Using  the  relations  V  x  A  =  B  and  32A  =  -  u  J 
in  the  expression  for  the  magnetic  energy  of  the  system 


-  / 


y  B' 


and  after  integrating  by  parts,  we  find 

e  =  -  UQ  Jd3r  J •  A  =  -  iryo  JrJ0A0drd2.  (5) 

Substituting  in  Eq.  (5)  the  Fourier  expansion  of  3  and  A, 
it  is  obtained 

N-l  N-l  N-l 

e  -  -wp0  /dr  r£  ),.£  a,  £  e  ~  ,  (6) 

k=0  1=0  m=0 

where  N  is  the  number  of  axial  cells  and  the  integral  over 
z  has  been  reduced  to  a  sum  over  m.  Using  the  identity 


2-rriLm 

2^  e  N  =  N,  for  L  =  0  and  L  =  N 
m=0  =  0,  otherwise 


Eq.  (6)  becomes 


N-l 

e  -  -  "“o  Jdr  rZ  ^  aN-k’ 

k=0 

Because  both  J  and  A  are  real  quantities 


a 


N-k 


* 
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and 

.  * 

jN-k  =  jk  * 


Substituting  these  two  relations  into  Eq.  (7) ,  the  magnetic 
energy  takes  the  following  form 
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-7TU 


0 
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k=l 


rdr , 


or 


£ 


N/2 


where 


e  =  -  2*UqAz  /  (jk  rdr. 
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m 


and 


e  =  -  ttu  dz  f  j  a  rdr. 
o  o  J  Jo  o 

The  quantities  may  be  identified  as  the  magnetic  energy 

associated  with  the  k  -  th  eigenmode  of  the  system.  Initially 
most  of  the  magnetic  energy  in  a  layer  will  be  in  eQ  because 
initially  Jg  and  Ag  do  not  vary  in  z .  In  unstable  systems 
the  dominate  non- zero  mode  corresponds  to  the  number  of  rings 
that  the  layer  has  formed. 

III.  THEORETICAL  PREDICTIONS 

a)  Stabilization  with  an  Axial  Velocity  Spread 
The  tearing  mode  stability  properties  of  a  strong  ion 
layer  immersed  in  a  background  plasma  have  been  studied 
recently  by  Uhm  and  Davidson.7  Using  a  hybrid  model,  in 
which  the  thin  ion  layer  is  treated  kinetically  and  the  back¬ 
ground  plasma  as  a  fluid,  they  have  investigated  perturbations 
with  frequencies  jwj<  w  .  ,  where  <£  .  is  the  ion  cyclotron  fre- 
quency  at  the  outer  edge  of  the  layer.  When  the  radial  thick¬ 
ness  2a  of  the  layer  is  much  smaller  than  its  mean  ’■adius  Rc, 
they  have  shown  that  perturbations  are  stabilized  when  the 
axial  velocity  spread  dvz  of  the  ions  in  the  rotating  layer 
exceeds  a  critical  velocity  spread  given  by 


i 


i 
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dV 


crit 


(R„u3  •  ) 

O  Cl 


-  [l  -  ‘W’]  Tifer  (x  *  357-) 
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where  Rc  is  the  radius  of  conducting  wall  surrounding  the 
layer  and  v  is  the  Budker  parameter. 


An  important  feature  of  Eq.  (8)  is  that  the  critical 
axial  velocity  spread  required  for  stabilization  is  independent 
of  the  background  plasma  and  therefore  should  be  valid  even  when 
the  plasma  density  goes  to  zero  as  in  the  computer  simulation, 
b)  Wall  Stabilization 

A  conducting  outer  wall  plays  an  important  role  in  the 
tearing  mode.  The  expression  for  hv  .  shows  that  when 

wi  J,  L 

Rq  =  Rc  the  system  is  stable  to  tearing  instability.  The 
behavior  of  a  system  in  which  all  parameters  except  the  wall 
radius  are  held  constant  is  much  more  complex  than  equation  (8) 
would  indicate  at  first  sight.  This  is  due  to  uici  being  a 
function  of  the  wall  radius.  Specifically 

u>ci  =  -  (\>+l)wco  [v(2  x2-  1)  -  l]'1  ,  (9) 

where  uco  is  the  ion  cyclotron  frequency  in  the  external 
magnetic  field  BQ  and  x  =  R0/&c* 

After  substitution  of  Eq.  (9)  into  Eq .  (8)  we.  obtain 

iv2  ..  *  V2F  (x)  (10) 


where 


fix,  =  -i-1— xi>..v - , 

[v ( 2x2 -1) -1] 2 


and 

V2  m{1  +  R°  “co  • 

V  is  a  function  of  v  but  is  independent  of  the  wall  radius. 

In  the  limit  of  a  thin  beam,  i.e.,  for  a  0,V  is  independent 

of  u.  F(x)  is  a  complicated  function  containing  a  singularity. 
Figure  1  shows  the  dependence  of  F(x)  on  RQ/RC  for  v  -  0.3, 

0.8  and  2.0.  Because  of  the  singularity  F(x)  is  artifically 
suppressed  in  the  plot. 

R  /R  =  0  corresponds  to  no  wall  or  the  limit  as  the 
o  c 

wall  is  removed  to  infinity.  The  other  extreme  RQ/RC  =  1 

corresponds  to  a  very  thin  beam  against  the  conducting  wall. 

The  lower  the  value  of  F(x)  the  less  the  thermal  spread 

required  to  obtain  a  stable  layer.  As  shown  in  Fig.  1  the 

function  F(x)  is  characterized  by  three  distinct  regions. 

For  v  <  1/3,  F(x)  is  monotonically  decreasing  as  RQ/RC 

increases.  Thus,  if  a  system  is  stable  with  no  conducting 

wall  (i.e.,  R_/R_  =  0)  it  is  stable  for  all  values  of  R  /R  . 

o  c  o  c 

For  1/3  <  v  <  1,  F  has  a  maximum  and  therefore  the  system 
requires  a  greater  axial  velocity  spread  to  be  stabilized 
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as  the  conducting  wall  moves  from  far  away  closer  to  the 
layer.  For  values  of  (RQ/R  )  greater  than  that  corresponding 
to  the  peak  of  F(x) ,  the  equilibrium  becomes  more  stable. 

For  v  >  1,  the  function  F(x)  has  a  singularity  at 
(RQ/Rc)s^n  =  (1  +  v)/2v  .  The  conducting  wall  has  a  destabi¬ 
lizing  effect  on  the  layer  when  its  radius  is  such  that  the 
condition  RQ/RC  <  (Rc/Rc^sin  satis^ied-  In  contrast,  when 
the  wall  is  located  very  close  to  the  outer  surface  of  the 
layer,  i.e. ,  where  RQ/RC  >  (R0/Rc) sin  the  stability  of  the 
system  is  improved.  Note  that  at  the  singularity  equilibrium 
exist  only  when  the  external  magnetic  field  is  equal  to  zero. 
IV.  RESULTS  OF  SIMULATIONS 

Simulations  using  6  to  24  thousand  macro-particles  show 
excellent  agreement  with  the  predictions  of  Eq.  (8)  even  when 
the  inequality  a  <<  RQ  is  not  satisfied, provided  that  Rq  is 
set  equal  to  |"(R2  +  R2)/^"]*5,  where  R  is  the  inner  and  R  is 
the  outer  radii  of  the  layer. 

The  noise  level  of  the  system  varies  as  (number  of 

-k 

particles)  .  Since  the  number  of  particles  in  the  simulation 
is  limited,  even  a  stable  system  for  which  the  linear  theory 
would  predict  no  growth,  the  axial  components  of  the  Fourier 
spectrum  will  be  non-zero  with  the  equilibrium  fluctuation 
level  determined  by  the  number  of  particles  used  in  the 
simulation. 
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Although  some  degree  of  randomness  is  provided  by  the 

initialization  process,  it  is  well  below  the  final  thermal 

levels  found  in  stable  equilibria.  For  this  reason  each 

plot  of  the  energy  in  the  eigenmodes  is  characterized  by  an 

initially  fast  rising  segment  which  is  terminated  by  reaching 

the  thermal  level,  usually  after  one  angular  rotation  period. 

Figure  2  shows  the  noise  levels  vs.  the  normalized  time 

tw  =  tw  .  (v+1)-1  from  two  simulations  with  6000  and  24000 

R  Cl 

macro-particles  stabilized  by  axial  temperature.  The  two 
simulations  are  identical  except  for  the  number  of  particles. 

The  noise  drops  by  about  a  factor  of  two  when  four  times  as 
many  particles  are  used. 

Figure  3  shows  details  from  the  simulation  of  6000  macro¬ 
particles.  In  this  run  the  Budker  parameter  is  2,  the  inner 
and  outer  radii  of  the  layer  are  20  and  26  centimeters  respec¬ 
tively,  the  wall  radius  is  32  centimeters,  the  applied  magnetic 
field  is  10  kilogauss  and  the  thermal  spread  is  1.2  times  greater 
than  the  critical  value.  The  periodic  temporal  variation  of  the 
individual  axial  components  corresponds  roughly  to  the  time 
required  for  a  particle  with  average  axial  velocity  to  move  a 
wavelength. 

If  the  thermal  spread  Av  is  less  than  Avcrj_t'  the  modes 

grow  with  time.  A  comparison  of  the  growth  rates  of  individual 

14 

modes  with  those  predicted  by  linear  theory  for  np  = 
a  discrepancy  of  about  a  factor  of  two,  with  the  code 


0  shows 


consistently  giving  lower  growth  rates.  Figure  4  shows  the 
time  development  of  an  unstable  system.  All  the  parameters 
are  identical  to  the  stable  case  except  that  the  initial 
axial  thermal  spread  is  a  quarter  of  the  value  required  for 
stability  i.e.,  Av  =  0.25  Av  . , .  The  measurement  of  growth 
rates  is  complicated  by  the  existence  of  the  ambient  noise. 

The  growth  rates  are  compared  to  the  theory  when  the  long  wave¬ 
length  axial  modes  have  established  a  measurable  growth  rate. 

As  the  layer  tears  the  axial  distribution  function 
changes  from  the  initially  square  configuration  to  a  more 
Maxwellian  shape.  Figure  5  shows  the  rather  dramatic  broadening 
that  has  occurred  by  the  time  the  growth  rates  are  measured. 

The  broader  profile  probably  accounts  for  the  discrepancy 
between  our  measurements  and  the  theory.  The  broadening  is 
indicative  of  the  bulk  heating  of  the  layer.  The  unstable 
spectrum  initially  grows  for  short  wavelengths  and  as  the 
heating  stabilizes  these  modes  the  evolution  continues  towards 
longer  wavelengths. 

The  correlation  between  the  energy  in  the  dominant 
eigenmode  and  the  number  of  rings  that  a  layer  has  broken  into 
is  illustrated  in  Fig.  6.  As  shown  in  Fig.  4  the  system  is 
dominated  at  t  =  6  ul1  by  the  second  eigenmode.  As  expected 

X\ 

the  system  has  formed  two  well  defined  particle  rings  at  this 
time.  In  the  final  stages  of  development  the  layer  is  dominated 


by  the  one  mode  and  consists  of  a  single  ring.  An  examination 
of  the  axial  velocity  distribution  reveals  that  in  this  fixed 
state  that  the  particles  have  a  spread  greater  than  the  criti¬ 
cal  spread  required  for  stabilization. 

To  verify  the  possible  destabilizing  influence  of  a 
conducting  wall,  simulations  of  two  nearly  identical  but  care¬ 
fully  selected  systems  were  made.  The  only  external  difference 
between  the  two  systems  was  the  radius  of  the  conducting  outer 
wall.  The  layer  extended  from  8.8cm  to  16.8cm,  was  immersed 
in  an  external  axial  magnetic  field  of  10  kilogauss  and  had  a 
Budker  parameter  of  2.  With  the  outer  wall  at  32cm  an  axial 
thermal,  spread  was  selected  to  give  stability.  This  same 
thermal  spread,  however,  was  not  sufficient  to  stabilize  the 
layer  if  the  wall  radius  was  reduced  to  21cm. 

The  two  cases  described  corresDond  to  R  /R  =0.4  and 

o  c 

0.6  for  a  Budker  parameter  greater  than  one  as  seen  in  Fig.  1. 

A  time  history  of  the  energy  in  the  axial  structure  of  the  two 
systems  is  shown  in  Fig.  7.  For  the  case  RQ/RC  =  0.4  (the 
more  distant  wall)  the  layer  developes  to  the  thermal  level 
and  does  not  grow.  For  R0/Rc  =  0.6  (the  close  wall)  the  layer 
exhibits  growth  in  the  axial  modes  well  above  the  thermal  level 
and  tears. 

V.  SUMMARY  AND  CONCLUSIONS 

These  computer  simulation  experiments  have  verified  the 
predictions  of  Uhm  and  Davidson's  theory  about  the  critical 
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axial  velocity  spread  required  to  stabilize  the  tearing  mode 
in  an  infinitely  long  ion  layer.  Although  the  theory  has 
developed  under  the  assumption  that  :he  thickness  of  the 
layer  is  much  smaller  than  its  average  radius,  the  results 
of  the  simulation  have  shown  that  the  predictions  of  the 
theory  are  even  valid  when  a  <<  R  ,  provided  that  Rq  is  set 
equal  to  the  geometric  mean  radius  of  the  layer. 

The  most  important  conclusion  of  the  present  work  is 
on  the  effect  a  conducting  tube  surrounding  the  layer  has 
on  its  stability.  For  a  very  long  time,  it  has  been  commonly 
accepted  that  a  conducting  waveguide  improves  the  tearing  mode 
stability  properties  of  the  layer.  However,  our  results  show 
that  when  v  >  1/3,  a  conducting  wall  may  have  destabilizing 
influence  on  the  layer. 
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Fig.  1  -  Dependence  of  critical  thermal  spread  on  wall  radius.  Solid  line  is 
for  a  Budker  parameter  less  than  1/3.  Dashed  line  is  for  a  Budker  parameter 
greater  than  1/3  and  less  than  1.  Chain-'dot  is  for  a  Budker  parameter  great¬ 
er  than  1. 
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Fig.  2  -  Noise  levels  In  identical  systems  with  different  numbers  of  macro¬ 
particles.  Upper  curve  corresponds  to  about  6,000  and  lower  curve  corresponds 
to  24,000  macro-particles. 
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Energy  in  axial  Fourier  spectrum  vs  time  in  a  system  stabilized  by 
axial  thermal  spread.  The  parameters  are  given  in  the  text. 
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